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Abstract 

Recent single-molecule pulling experiments have shown how it is possible to ma- 
nipulate RNA molecules using optical tweezers force microscopy. We investigate a 
minimal model for the experimental setup which includes a RNA molecule connected 
to two polymers (handles) and a bead, trapped in the optical potential, attached to one 
of the handles. Initially, we focus on small single-domain RNA molecules which unfold 
in a cooperative way. The model qualitatively reproduces the experimental results and 
allow us to investigate the influence of the bead and handles on the unfolding reaction. 
A main ingredient of our model is to consider the appropriate statistical ensemble and 
the corresponding thermodynamic potential describing thermal fluctuations in the 
system. We then investigate several questions relevant to extract thermodynamic in- 
formation from the experimental data . Next, we study the kinetics using a dynamical 
model. Finally, we address the more general problem of a multidomain RNA molecule 
with Af g 2+ -tertiary contacts that unfolds in a sequential way and propose techniques 
to analyze the breakage force data in order to obtain the reliable kinetics parameters 
that characterize each domain. 

1 Introduction 

The RNA molecule plays a central role in molecular biology showing an enzymatic 
function during the translation and splicing processes E] • Experiments based on 
the manipulation of single-biomolecules, such as laser tweezers force microscopy, allow 
scientists to investigate their mechanical properties. These give information about the 
structure, stability and the interactions involved in the formation of such structures 
pU El El E| ■ In these experiments mechanical force is applied to the ends of 
a RNA molecule. The molecule is then pulled |l()l ITT] until a value of the force is 
reached such that the molecule unfolds. If the pulling process is reversed then the 
molecule refolds again. In these experiments the force exerted upon the system is 
recorded as a function of the end-to-end distance giving the so-called force-extension 
curve (FEC). The nature of this unfolding-refolding process is stochastic and therefore 
the values of the force at which the molecule unfolds-refolds change from experiment 
to experiment. Sometimes, as in the case of presence of Mg 2+ -tertiary contacts, it 
is not possible to pull the molecule in quasi-static conditions because the relaxation 
time is too large for the experimental possibilities which are largely limited due to 
the presence of strong drift effects in the machine. Therefore, during the pulling 
process, the molecule is driven to a non-equilibrium state which is characterized by 
strong irreversibility effects. The study of this pulling process might be useful to 
understand many biological processes where biomolecules are unfolded under locally 
applied force, for example when the mRNA goes through the ribosome during the 
translation process. 

To manipulate a RNA molecule some synthesized polymers typically several hun- 
dred nanometers long (called handles) have to be chemically linked to the extremes 
of the RNA molecule. A polysterene bead is then chemically attached to the end of 



one of these handles and used to measure the force by reading its position inside the 
optical trap. These additional elements (bead and handles) are an inseparable part of 
any pulling experiment and they have an influence on the unfolding process. To char- 
acterize the thermal behavior of the pulled global system (bead, handles plus RNA 
molecule) it is important to identify the proper control parameter. This is an essen- 
tial step towards the modclization of the experiment and has several consequences. 
For instance, the force acting on the extremes of the RNA molecule cannot be ex- 
ternally controlled but fluctuates and its mean value depends in a non-linear way on 
the value of the control parameter. The control parameter determines the relevant 
thermodynamic potential defining the equilibrium state of the global system as well 
as the magnitude of the fluctuations around that state. A proper inclusion of these 
parts is necessary to accurately interpret the experimental data. Another important 
point of the work is the model for the RNA molecule. We consider the RNA molecule 
to be composed by different domains, each one showing cooperative unfolding. Each 
domain is then modeled as a two-states system: the unfolded state (UF) and the 
folded one (F), which are separated by a kinetic barrier. A main effort throughout 
this paper is to present in the most clear way the appropriate theoretical frame to 
understand pulling experiments leaving aside further additional complications, never- 
theless important, such as a detailed response of the optical tweezers machine or the 
microscopic structure of the RNA molecule. 

The goal of this paper is twofold: (i) we show how to build a minimal model aiming 
to reproduce the experimental setup including all the aforementioned elements (bead, 
handles and the RNA molecule) and quantitatively reproducing various experimental 
results; (ii) we show how to analyze experimental data extracted from both quasi- 
static and out-of-equilibrium pulling experiments in order to obtain thermodynamic 
and kinetic information about the unfolding reaction. 

The paper is divided into three main parts. In the first part of the paper (Sec- 
tions EGGJ we describe the model for the experimental setup (Sec. EJ and introduce 
the ensemble that is relevant to model the pulling experiment fScc. 12. lj) . In Sec. El 
we describe the two-states model convenient to reproduce the cooperative unfolding 
of the RNA molecule and in Sec. 0] we describe the models used for the bead and han- 
dles. In the second part of the paper (Sections I5I6|I we analyze the unfolding-refolding 
behavior of a cooperative two-states RNA molecule in a pulling experiment for both 
equilibrium and non-equilibrium regimes. For the equilibrium regime, we compute 
the partition function in the ensemble that is experimentally relevant, and derive an 
expression for the quasi-static work exerted upon the system as the molecules unfolds. 
This expression relates the work measured in a quasi-static pulling process to the dif- 
ference of free energy between the F and UF states at zero force AG . We analyze 
in detail the different thermodynamic contributions to the total work, the influence 
of the parameters describing bead and handles on the FEC, and obtain expressions 
for the force at the midpoint of the transition. For the non-equilibrium behavior we 
investigate in detail the fraction of molecules that unfold (refold) more than once 
during the unfolding (refolding) path, which is a quantity amenable to experimental 
checks. We find that this fraction is related to the mean dissipated work exerted upon 
the system, which gives us a way to extract the reversible work in non-equilibrium 
processes just by measuring the total work. We also identify an interesting symmetry 
property relating these fractions for the forward and reverse processes. To endorse 
most of our theoretical results we also consider a simulation of a pulling experiment 
that allow us to obtain the characteristic FEC, cither in a situation where the tran- 
sition occurs in equilibrium or in a situation where it does not. In the third part of 
the paper (Sec.[7|), we address the unfolding behavior of complex RNA molecules with 
more than one folded-domain and in the presence of Mg 2+ -dependent barriers. In 
this case, the refolding is not observed at the experimental conditions, and the dis- 
tribution of the breakage force is a first order Markov process Q3| ■ We focus our 
attention in the specific case of RNA molecules where domains unfold in a sequential 
fashion according to a reproducible path. This unfolding mechanism is generally a 
consequence of the topological connectivity of the different parts of the molecule and 
of the blockade of the force induced by the most external tertiary contacts on the 
interior domains. We model the molecule as a series of domains, each represented 
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by a two-states system, and we compute the distribution of breakage force for each 
domain. We propose several methods of analyzing the breakage force data in order 
to achieve reliable values for the height and position of the barrier of each domain. 
In Sec. ISJwe present the conclusions. Five appendixes are devoted to describe some 
analytical calculations. 



2 Model for the experimental setup 

We consider a minimal model in order to reproduce the experimental setup of a pulling 
experiment carried out using laser tweezers force microscopy |101 114j . The model 
(Fig. Q is composed by a small RNA molecule connected to two polymers called 
handles 1 which are used to attach the small RNA molecule to two beads at each end. 
One bead is confined in the optical trap generated by the laser beams, the other is 
held fixed to the tip of a micropipette by air suction. 




Figure 1: Schematic picture of the model for the experimental setup in a RNA pulling experiment 
as described in the text. We show the configurational variables of the system xj,,x r , x% x , Xh 2 which 
are the projections of the end-to-end distance of each element along the reaction coordinate axis 
(i.e. the axis along which the force is applied). The potential Vb(xb) is well described by an 
harmonic potential of one-dimensional spring with rest position at Xb = 0. 



The whole system consists of a chain of connected elements. Starting from the 
left side of the chain there is a bead (i?i) of radius i?bead that is trapped in the laser 
tweezers potential, Vb(x) 2 . We use the position of this bead B\ to read the force 
acting on the system in the same way as the needle of a 'manometer' is used to read 
the pressure exerted by a gas on the walls of a container 3 . The second element is 
a handle (handle 1) with one end specifically attached to the bead B\ and the other 
end attached to the RNA molecule at the point A. The second handle (handle 2) has 
one end specifically attached to the RNA molecule at the point C. The other end is 
specifically attached to the bead B2, fixed to the tip of a micropipette. The molecule is 
pulled by moving the micropipette along the x direction . The configurational variables 
of this simplified system are taken as the projections of the end-to-end distances of 
each element along the force axis: Xh ± = B\A — Rbe&d, Xh 2 = CB2 — i?bead for the 
distances of the handles, x r — AC for the RNA end-to-end distance and x& for the 

1 These are hybrids of DNA and RNA rather than single stranded DNA or RNA polymers in order to 
avoid the formation of secondary structures 

2 Although the trap potential should be defined in the three-dimensional space Vb(x) we will consider 
Vb(x) as the potential of mean-force projected along the reaction coordinate axis. This approximation is 
very accurate as fluctuations along the y, z directions are assumed to hardly affect the unfolding behavior 
of the RNA molecule. 

3 This is not the way the force is usually measured in dual beam optical tweezers where two photodetec- 
tors located at opposite sides of the chamber are used to collect the total amount of deflected light which 
is then converted into force after calibration of the machine. 
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position of the bead B\ in the trap. 
Xb of the bead B\: 

!-- 



The force / is measured by reading the position 
dV b (y) 



dtj 



(1) 



We define the subsystem S as that composed by the two handles and the small RNA 
molecule. The end-to-end distance for the subsystem S is then given by x = xi n + 
Xh 2 + x r (Fig- GJ- The total distance between the center of the trap and the tip of 
the micropipette is given by Xt + -Rboad = Xb + x + i?boad- Pulling experiments give 
FECs, f(x), corresponding to the force exerted on the chain (JJJ measured through 
the position of the bead B\ as a function of the end-to-end distance of subsystem S. 



2.1 Ensembles 



Experimentally it is possible to consider two different ensembles depending on which 
variable is used as the externally imposed non- fluctuating parameter 4 . 

• Mixed ensemble: The total distance between the center of the trap and the tip 
of the micropipette is held fixed, hence Xt is the externally controlled parameter. 
In this ensemble there are fluctuations in x and / given by |16l I17j , 



{Sx 2 } = 



k R T 



k x (X T ) + k b (X T ) 
with k x (Xx) 



(Sf) = 



k B Tkl 



d(f) 



d{x) 



k x {X T ) + k b {X T ) ' 
d(f) 



h(X 7 



d(x b ) 



(2) 



where (...) stands for thermal average, ks is the Boltzmann constant, T is the 
temperature of the bath, kb{Xx) is the stiffness of the optical trap and k x (Xx) is 
the effective rigidity corresponding to the subsystem S. The latter is determined 
by the serial compliance 



k x (X T ) 



1 



1 



+ 



1 



k hl (X T ) k h2 (X T ) k r {X T ) 



(3) 



where kh t (i = 1,2) and k r are the rigidities of the handles 1, 2 and the RNA 
respectively. These rigidities are Xt dependent and so are the fluctuations J5J. 

Force ensemble: In this case a piezo actuator controls the force (and therefore 
the position of the bead Bi). In this ensemble Xt and x are fluctuating variables, 
(SX T ) — (5x 2 ) — ksT /k x {f), where k x (f) is the stiffness of the subsystem S 



when the force is held fixed, k x (f) 



df 



Most of the theoretical work for the denaturation of RNA in pulling experiments 
considers the force ensemble. However, it is experimentally very difficult to work in 
the force ensemble where either the force or the variable Xb must be controlled. Indeed, 
for Xt to fluctuate the center of the trap must also fluctuate to compensate for the 
fluctuations in the force. It is difficult to imagine how to experimentally implement 
such ensemble. Therefore the most natural ensemble is that where Xt is constant. 
Indeed this is the ensemble most relevant for the experiments and therefore we will 
work in the mixed ensemble throughout this article. 



3 Two-states model for a single RNA domain under 
the effect of an external mechanical force 

The unfolding of some biomolecules under the effect of a mechanical force is a highly 
cooperative process that can be qualitatively described by a two-states model. The 
two-states model has a long tradition in physics and has been applied previously by 
several authors in order to explain the unfolding behavior of single domains of proteins 

4 The existence of an external non-fluctuating parameter is required to have a well defined equilibrium 
state. 
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and RNA hairpins [H] EE1 EH EDI EH E21 Recently, it has been shown how such a 
simple phenomenological description, with Kramer transitions-rates, does not fully 
reproduce the kinetics observed in pulling experiments of the protein Titin, and more 
realistic descriptions have been proposed [231 • 

Let us consider an RNA molecule isolated from the rest of the system in equilibrium 
at constant temperature, pressure and zero force. In the simplest description both 
states (hereafter denoted by UF -unfolded- and F -folded-) are characterized by their 
Gibbs free energy Gy F and G F respectively and the RNA molecule occupies each 
state with a probability given by the Boltzmann distribution. In a more refined 
description the molecule can also occupy intermediate configurations depending on 
the number n of the first-opened, or denaturated, bases ( Fig. [3 (a)) [2I]. The F and 
UF states correspond then to the RNA configuration with n — and n — N bases 
opened, where N is the total number of pair of bases of the molecule. The free energy 
landscape is described by a function G°(n) which characterizes the probability of a 
hairpin having the first n bases opened (Fig. [21 (b)). This description excludes the 
existence of other breathing intermediate configurations that might be relevant for the 
unfolding reaction [25| . 




tl = 8 (number of opened bases) 



(a) I n I 

(b) 

Figure 2: Schematic representation of (a) a RNA hairpin with n = 8 bases opened, (b) the free 
energy landscape for a single RNA hairpin at zero force as a function of the number of denaturated 
bases n at T < T me iti n g (melting temperature of the RNA) and normal ionic conditions. In this 
situation the stable state is the folded one with n = Q. 



When an externally controlled force / is applied to the ends of the RNA molecule 
the adequate thermodynamic potential to consider is the Legendre transform of the 
Gibbs free energy G'(n) = G°(n) — fx r (n) (2SJ. The free energy landscape G' is then 
tilted along the reaction coordinate x r , which is the projection of the end-to-end dis- 
tance of the molecule in the axis force and explicitly depends on the number of opened 
bases n. Since we work in the ensemble where neither / nor x r are control parameters 
the non-fluctuating parameter Xt determines the adequate thermodynamic potential 
Gx T ■ The free energy Gx T of the system shown in Fig.^is a potential of a mean force 
that characterizes the equilibrium state of the whole system, including the handles, 
the bead and the RNA molecule, at fixed value of Xt- In order to build the model is 
useful to represent the free energy Gx T i as a function of the end-to-end distance of 
the subsystem S, as shown in Fig. [31 (a). This picture tells us about the probability 
Px T (%) of finding the subsystem S at a given value of its end-to-end distance x for a 
fixed value of Xt, px T ( x ) x exp(— f3Gx T ( x ))> where f3 = 1/ksT. 
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(a) 

Figure 3: (a) Schematic representation of the free energy landscape Gx T ( x ) f° r the whole system 
at T < T mc itin g , for Xt < X T (where X T is the value of Xt where both states F and UF are 
equiprobable) and normal ionic conditions. In this situation the stable state is still the folded one. 
We also show all the parameters characterizing the two states model, (b) Schematic representation 
of the relevant configurations in the F and UF states along the reaction coordinate x r . We consider 
the F state to be characterized by a single configuration x r — and the UF state by a continuous 
set of values of x r . We use the label a = for the F state and a = 1 for the UF state. 

The free energy landscape Gx T (x) shows two pronounced minima corresponding 
to the F and UF states (Fig. |3J). The discrete variable a stands for the state of the 
domain: the value a — denotes the F state and a = 1 the UF state. The relative 
thermodynamic stability of these states depends on the difference of free energy be- 
tween them, AG(Xt)- Moreover we will consider the existence of an intermediate or 
transition state along the reaction path from the F to the UF states and vice versa. 
This transition state is the intermediate RNA state with highest free energy connect- 
ing the F and the UF sate along the reaction path. It may correspond to a RNA 
configuration where the first n — n* bases are opened 5 . In the simplest scenario the 
intermediate state can be assumed to have a very short lifetime, its main effect is to 
hinder transitions from the F to the UF state and back. In this scenario the transition 
state can be represented by an activation barrier and this is the model we will adopt 
throughout the paper. The F and UF states are separated by a barrier of height 
B(Xt) measured relative to the F state. The barrier is located at a distance x\{Xt) 
from the F state and X2(Xt) from the UF state. The distance between the two states 
is i m (lr) = x\{Xt) + Xz{Xt)- Because the rigidity of the RNA molecule in the F 
state is very large we can assume this state to be characterized by a single configura- 
tion corresponding to the value x r = of the reaction coordinate; fluctuations around 
this configuration cost so much energy that they are highly improbable. The RNA 
in the UF state has a finite rigidity, hence it is represented by a set of configurations 
within a continuous range of values of x r (Fig. [31 (b)) ■ 

4 Modeling the different parts of the setup 

In this section we specify how we model the different elements of the system: the bead 
trapped in the optical tweezers potential, the two handles, and the RNA. 

We stress that the shape of the free energy landscape depends on Xt as well as the location of the 
barrier corresponding to the intermediate state. Therefore the value of n* that characterizes the transition 
state is also Xt dependent. There are experimental limitations to follow the folding and unfolding of the 
molecule (hopping) given by the operational range of frequencies of the instrument used. For instance, 
in this operational range was 0.05 — 20Hz, meaning that hopping events out of this range were not 
observable. Moreover, in pulling processes the folding-unfolding reaction only occurs in a narrow range 
of values of the control parameter Xt around X T , otherwise the folding- unfolding relaxation time is too 
large. Hence the study of the folding-unfolding kinetics is restricted to the regime Xt ~ X T and to the 
operational range of frequencies. Therefore, as the transition state n* is only relevant for the study of the 
kinetics, we assume n* independent of Xt, n* = n*(X T )- 
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4.1 Model for optical tweezers: a bead matched to a spring 



Typical experimental values for the trap stiffness and the diameter of the beads are 
kb w 0.15 — 0.05pN/nm and -Rbcad ~ 1 — 3/im respectively. We consider that the bead 
follows a Langevin dynamics of an overdamped particle (i.e. without inertial term) 6 : 

^ = F R (x b ) + m , (4) 

where 7 is the friction coefficient and Fr is the resultant force applied to the bead. 
The stochastic term is a white noise with mean value (£(£)) = and variance 
= 2ksT^S(t — t'). The force Fr has two contributions: one coming from 
the optical trap potential, /, given by lfl]l. and the other from the subsystem S, 
f x 7 . Therefore Fr = —f + f x (xb). The experimental results |10l show a linear 
dependence of the force / on x along the transition (rip) where the variable x refers to 
the subsystem S (Fig.^J, hence we conclude that the optical trap V b is well modeled 
by an harmonic potential of stiffness kb- We can then express the force / measured 
through the optical tweezers as / = fc&Xf, = kb(Xx — x), where we have used In 
equilibrium the average force acting upon the bead is zero (Fr), hence (f x ) = kb(xb}- 
However Xb fluctuates and so both instantaneous forces / and f x are not identical. 
Doing an expansion around the equilibrium position of the bead, x eq 8 , we obtain: 

7"^- = -kR(x b - x eq ) + , (5) 

where kR is the effective spring constant applied to the bead, kR = k x + kb, with k x 
given by J2J. The relaxation time of the system (i.e the typical time during which the 
position of the bead decorrelates), Tb, is given by Tb — j/kR. Applying the Stokes's 
law for the friction coefficient in water we obtain: 7 = QitRr) w 10 _5 pNs/nm. The 
stiffness of the handles and the RNA are force dependent. Near the F-UF transition, 
typically these stiffness values are, at least, one order of magnitude bigger than kb- 
Taking kb = O.lpN/nm, and kh x , kh 2 , k a > lpN/nm we get Tb < 10~ 5 s 9 . By collecting 
data at frequencies smaller than 10 5 Hz we can guarantee that we will not have effects 
due to the bead's overdamping, hence the bead relaxes quickly to its new equilibrium 
position at each step. This ensures that assuming an instantaneous relaxation of the 
bead position is enough to capture its overdamped dynamics. 



4.2 Polymer model for handles and single-stranded RNA 

The handles and the single-stranded RNA (ssRNA), the UF state of the RNA, are 
polymers that typically measure d ~ 1 — 3nm in diameter and L w 20 — 500nm in 
length. As the bead has a much bigger size than the polymers the friction coefficient 
(and therefore the relaxation time) for the polymers is much smaller. This allows us 
(as we do for the bead) to only consider an instantaneous relaxation for the handles 
and the ssRNA. To describe the equilibrium behavior of the handles and the ssRNA 
under the effect of an external force we use the worm-like-chain (WLC) model |27j . 
This is described by a Hamiltonian that includes a local bending term as well as the 
potential energy of the polymer in the presence of the pulling force. Parameterizing 
the polymer with the arc length s, the energy function can be written as: 



[ L ° \k B TP fdt(s)\ 2 „ 



K 1 > 1 : i 1 1 i \ •< \ \ ~ 

ds , (6) 



6 In Q| we are neglecting the drag force felt by the bead (equal to —71;) as the chamber is moved (and 
the water dragged relative to the lab frame) at a certain pulling speed v = —y^. For the range of pulling 
speeds used in the experiments this contribution is negligible, of the order of O.lpN. 

7 This is also the force exerted upon the subsystem 5* for a given value of x = Xt — Xb. 

8 We expand f x and / around x eq keeping only the first term in (xb — x eq ), i.e f x ~ (f x ) + k x {x eq — xt), 
with k x given by ©, and / ~ (/} + kt(xb — x eq ). This approximation is valid in the regime where 
the force fluctuations are not big. Using that at equilibrium {/} = (f x ) we obtain Fr = — / + f x ~ 
-{k x + k b )(x b - x eq ). 

9 In absence of handles kn = k b and r b = 10~ 3 s. Therefore 10 -3 s is the slowest relaxation time for the 
bead corresponding to the regime where the handles have practically no rigidity, k x ~ 0, a situation only 
encountered at small forces (below lpN approximately). 
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where L a is the contour length of the polymer, t(s) is the unit tangent vector along 
s, 9(s) is the angle between t(s) and the force axis, and P is the persistence length 
defined as the typical distance over which ^correlations decay to zero: (t(s)t(s')) w 

a-a'l 

e p . The persistence length of a polymer depends on the ionic conditions |28j . 
and typical values are 50nm for double-stranded DNA (dsDNA) and lnm for ssRNA. 
The thermodynamic properties of this model cannot be exactly computed, yet there 
are useful extrapolation formulas. Bustamante et al. |29| have proposed a simple 
expression for the force as a function of mean end-to-end distance of the polymer x, 



f 



k B T\ 1 



P 



4(1 - x/L ) 



l/4 + x/L . (7) 



Eq.J7J gives the exact solution as x approaches either zero or L a and is accurate at 
least up to 90% in between. Bouchiat et al. PU] have given an expression with an 
accuracy of 99% by adding to J7J) a polynomial of seventh order. The WLC model 
works well only at low forces, in the so called entropic regime where the molecule 
behaves as an entropic spring. At high forces there is an enthalpic correction due 
to the fact that the bonds are stretched and the contour length L a increases. To 
incorporate this effect it is then enough to replace x/L Q by x/L — f / E y in (JJJ, where 
E y is the Young modulus of the polymer (typical values are E y rj 500 — 1500pN for 
DNA-RNA molecules). 

Throughout Sees. and El we analyze the unfolding dynamics and thermodynamics 
of a single hairpin of RNA aiming to reproduce the results obtained from a pulling 
experiment for the hairpin P5ab in lOmM Mg 2+ JUj. In Sec. we use the partition 
function analysis to individuate the different thermodynamic contributions to the total 
free energy or reversible work done upon the system. Next in Sec. [B] we do numerical 
simulations of a pulling experiment. 



5 Thermodynamic analysis 

In this section we use the tools of statistical mechanics to analyze the thermodynamics 
of the system represented in Fig. ^ Most of the analytical development is done in 
appendix El an d in Sec. 15.11 we give the main results. In Sec. 15.21 we show how to 
get the force-extension curve (FEC) , the value of the force at the midpoint of the F- 
UF transition F c , and the different contributions to the total reversible work coming 
from the different elements of the system. In Sec. 15.31 we derive an expression that 
relates the reversible work exerted upon the subsystem S across the transition with 
the difference of free energy between the F and UF states at zero force, AG . As this is 
an experimentally measurable quantity this procedure provides a way to estimate the 
unfolding free energy of the molecule, a quantity biologically relevant as it determines 
the fate of biochemical reactions. Finally in Sec. 15.41 we show how to apply these 
relations to a specific example. 

5.1 Definitions 

In equilibrium the observables x a and the conjugated forces f a with a = h\,h2,r,b 
(referring to the different elements, handle 1 and 2, RNA and bead respectively) 
fluctuate. However, the thermodynamic free energy is only a function of the mean 
values of these observables that we denote by (x a ), (/ Q ). A representation of (f a ) 
versus (x a ) gives what we call the thermodynamic force extension curve (TFEC) for 
the element a in the mixed ensemble. If a refers to the whole subsystem S then the 
TFEC corresponds to the usual force-extension recorded in RNA pulling experiments, 
assuming that the pulling process is carried out reversibly. We can also define the 
restricted average (0),j(Xt) as the mean value of the observable O when the RNA 
molecule is in the state a for a fixed total end-to-end distance Xt- From now on, all the 
dependencies of the observables on the variable Xt will not be explicitly written, hence 
(0) a (X T ) = {0) a . In appendix El we derive an expression for the partition function 
Z(Xt), corresponding to the system schematically represented in Fig. ^ Applying 
the saddle point technique, and separating the contributions that come from the F 



8 



(<t = 0) and UF (a — 1) states we get: 



Z(X T ) = Z (X T ) + Z 1 {X T ) , 



where 



Z {X T ) « e X p[-(3(w hl ({x hl ) ) + W h2 ({x h2 ) ) + V b ((x b ) ) 



(8) 



(9) 



Zi(X r ) « exp[-)9(w hl ((a: /ll )i) + W h2 ((^ 2 )x) + H((a: 6 >i) + AG + W r ((a; r )i))] • (10) 

Here V& represents the optical trap potential and AG is the free energy difference 
between the F and the UF states at zero force. The function W a (x) corresponds to 
the reversible work performed by adiabatically stretching the element a from x a = 
to x a — x and it is given by 

W a (x a ) = / dxf a (x), with a = hi,h 2) r , (11) 
Jo 

where f a (x) is the TFEC for the element a (see appendix . We can define the 
probabilities for the RNA molecule of being in the F and the UF states by po and p\ 
respectively, 

Pa{X T ) = , with a = 0, 1 . (12) 

The thermodynamic value of any observable O can be expressed in terms of these 
probabilities, 

(O)= P0 (O) + Pl (O) 1 . (13) 

At the midpoint of the transition both states are equally probable, 

Po(X^) = Pl {X$) or Z (X^) = Z X {X%) , (14) 

where these functions have been defined in (|9I10() and (|12|) . Hence the midpoint of 
the transition in the mixed-ensemble is defined by the value of the control parameter 
Xf, that verifies 1(1*1)1. 



5.2 Computation of the transition force F c , the TFEC and the 
different contributions to the reversible work. 

The force at the transition, F c , is computed as the mean value of the force at X^ 
given by (|14fl . To reproduce the experimental results obtained for the P5ab RNA 
molecule in lOmM Mg 2+ JJO] we use the parameters given by Tables Hand H getting 
F c = 15.2pN 10 . This value is close to the one reported from the experiments F c c xp = 
14.5 ± lpN JHj- We also verify that the force at the transition F c is quite stable with 
respect to changes in the parameters of the problem used to model the handles and 
the bead trapped in the optical potential, such as the persistence and contour length 
of the handles, the spring constant and the bead radius. However, as the value of 
F c is highly influenced by the characteristics of the RNA molecule, we conclude that 
the dependence of the value of F c with the system is basically through the quantities 
AG , L r and P r . 



^r[pNnm] 


kb [pN / nm] 


Phi = P h2 [ nm ] 


Lh! = L h2 [nm] 




4.14 


0.1 


10 


160 


1000 



Table 1: Summary table of the parameter values used to model the handles and the bead in the 
optical trap. We use the value for the Young modulus corresponding to a dsDNA. The value for 
the other parameters have been taken from |10j. 

10 In our model we are considering the F state to be localized at x r = 0. However also the folded RNA has 
an end-to-end distance d r (the diameter of the hairpin) that tends to be aligned along the force axis when 
the force increases. Modeling the F state as a rigid segment of length d r = 2nm one obtains F c = 15.9pN 
which is a bit farther from the experimental value. 
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P r [nm] 


L r [nm] 




AG"[k B T] 


iV (number pair bases) 


1 


28.9 


800 


59 


22 



Table 2: Summary table of the parameter values used to model the RNA molecule. We use the 
value for the Young modulus corresponding to a ssDNA. The value for the other parameters have 
been taken from |1U| . 

Another interesting magnitude to measure is the reversible work W T ev done upon 
the system when pulling from an initial value of Xt = X T up to a final value of Xt- 
This work is given by 

W T CV (X T ) = G Xt - G x a = AG Xt .with 
Gx T = -k B T\n{Z{X T )) = -k B T\n{Z {X T ) + Z X {X T )) , (15) 

where we used |jSJ. The total reversible work in 1|15|) defines the change in the free en- 
ergy of the system. On the other hand the reversible work exerted upon the whole sys- 
tem is equal to the sum of reversible work exerted on each element Wjf v , W b ICV , W£ cv 
(handles 1 and 2, bead and RNA molecule) by changing the total end-to-end distance 
from the initial to the final value of Xt- 

W T CV (X T ) = Wl cv {X T ) + W h eY {X T ) + W* CV (X T ) , where (16) 
Wr(X T ) = (AV b )= Po (AV b ) + Pl (AV b ) 1 , (17) 

2 

wr(X T ) = (W h ) =J2[po(W hi ) +pi(W hi ) 1 ] , (18) 

i=l 

Wr(X T ) = (W r ) = Pi{(W r ) 1 + AG") , (19) 

where we used l|13|l . The functions AV b , Wh and W r correspond to the change in the 
potential energy of the bead in the optical trap and the work exerted upon both handles 
and the RNA molecule by moving the total end-to-end distance from the initial to the 
final value of Xt respectively. In Fig. 0] we show the different contributions to the 
total work TT™ V , W^ cv and W^ cv as a function of Xt as derived from the numerical 
computation of Z(Xt)- We also show the work WJP V exerted upon the whole system. 
Both computations H15I16J) overlap in a single curve as expected. Finally in Fig. 03 
(a) we represent the TFEC for the subsystem S, (/} n versus (x). This is obtained by 
numerical computation of the partition function using the relation, 

dG XT (X T ) d^ZjXr) i i y I \\ (irt 

We also present the results obtained by averaging 1000 different trajectories in a 
simulation of a pulling experiment as described later in Sec. and both curves show 
very good agreement. In Fig. [S] (b) we plot the mean force (/) as a function of the 
control parameter Xt- 



11 At a given state of the system (determined by a given value of Xt) all the forces f a are equal in 
average. Therefore the value of (/) referring to the mean force acting on the bead B\ (Fig. ^ coincides 
with the force (f a ) acting on each element a as well as on the subsystem S. 



10 



800 




Figure 4: Different contributions to the reversible work obtained from the partition function 
analysis: W T CV , W™ v , W b ICV and W* ev as a function of Xt- Note that the smallest contribution to 
the total work comes from the RNA molecule. 
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250 300 
<x> [nm] 

(a) 



350 




X T [nm] 

(b) 



Figure 5: (a) The continuous line corresponds to the results obtained from the numerical computa- 
tion of the TFEC. It is also shown the TFEC obtained by averaging over 1000 different trajectories 
in a simulation of a pulling experiment as explained in Sec. The pulling is carried out at an ap- 
proximate loading rate of 0.5pN/s (see footnote ll2|) . slow enough to generate a quasi-static process. 
One can observe that both curves agree, (b) Mean force (/) as a function of the control parameter 
Xt- Note that there is not an abrupt vertical drop of the force at X T . This is consequence of the 
narrow, yet observable, region of coexistence around the midpoint of the transition. 



5.3 Reversible work across the transition 

The quasi-static work Wf ip exerted upon the subsystem S across the transition is the 
area under the TFEC (Fig.©, from (x) = (x c ) to (x) = (x c ) u where the 

super-index c indicates that the system is at the midpoint of the transition where 
X T = X T JUJl, 

W T % = I*** h dy(f)(y) = V b (X T - (x c )t) - V b (X T - (x c ) ) . (21) 
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<f> 














W c 




<x> 



Figure 6: The shadow area under the TFEC along the transition corresponds to the quasi-static 
work W r 1 p (schematic representation). 

At the midpoint of the transition both states are equally populated and 114J1 holds. 
Therefore identifying © and we can write 1)21(1 as: 

W T % = AG + W? + AWH , (22) 

where the functions with a super-index c are evaluated at the mean value of their 
variables at the critical extension The W r is the loss of entropy of the RNA 

molecule along the transition due to the stretching and is given by (fill) , and the AW/, 
is the change of free energy of the handles across the transition computed as: 

AW h ^W hl {{x hl ) 1 ) + W h2 ((x h2 ) 1 )~W hl {{x hl ) )-W h2 {{x h2 ) ) . (23) 

Eq. (1221) tells us that the quasi-static work Wf ip coincides with the change of free 
energy of the different elements that form the subsystem S across the transition. This 
is experimentally measurable as the area under the rip observed in the TFEC 
corresponding to the F-UF transition (Fig. [fJJl . Therefore we can use (|22H to estimate 
AG from the TFEC, as explained in the next section. 

5.4 Estimate of AG from the TFEC 

In Fig.0we show two TFECs obtained from the partition function analysis correspond- 
ing to two systems with different kb but with the same handles and RNA molecule 
with parameters given in Tables ^andl^l respectively. We use (|22|l in order to extract 
the value of AG by computing W£ lp as the area under the rip in the TFEC (Fig. |BJ . 




<x> [nm] 



Figure 7: TFEC corresponding to two systems with handles and RNA characterized by the 
parameters given in Tables^[21and with an optical trap stiffness fct> = O.lpN/nm and fcf, = lpN/nm 
respectively. Note that the slope at the transition (rip) is proportional to — fcf,. 
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As expected for an harmonic trap l|2()(l . the TFEC in Fig. shows an slope at 
the transition (rip) proportional to — kb- To obtain the different contributions to 1)22(1 
using the WLC model [3D1 we nrs t estimate (x r )i and as the RNA and the 

handle i extension at the force value after the rip, respectively. In an analogous 
way we estimate the (27^)0 as the handle i extension at the force before the rip (/)o- 
Then we compute W£ and AW£ given by fTTJl and respectively using the WLC 
model [HOj for the TFEC f a (x) for the element a with a — hi,ti2,r (handles 1 and 
2 and RNA respectively). Finally, we compute the area under the TFEC across the 
transition (rip) in order to obtain and use (|2~2")l to extract AG . In Table we 
show the results obtained. 



fcb[pN/nm] 


W r c [k B T] 


AW£[k B T\ 


W r %[k B T] 


AG U [k B T] 


0.1 


20 


-8.5 


70.5 


59 


1 


17 


-41 


35 


59 



Table 3: Different contributions to the free energy change across the transition. As expected the 
value of AG is independent of the other parameters of the system. 

Note that the contribution AW 7 ^ is negative because when the RNA molecule opens 
the force relaxes and the handles contract, hence the free energy of the handles across 
the transition decreases. Neglecting the contribution that comes from the handles 
across the transition is a typical approximation often applied to experimental results. 
However we note here that this is not always possible as this contribution can be 
large. In the previous example, even in the case of small kb, we would loose 8fcsT in 
the balance equation (|22ll . The best condition to apply this approximation is to use 
handles characterized by a small ratio Lh/Ph as compared to the corresponding value 
for the RNA molecule {Lh/Ph « L r /P r ) and a potential well with stiffness as small 
as possible (i.e, small kb). In Fig. we show, for a small value of kb (kb = O.lpN/nm), 
how the different contributions to (1221 change when considering systems with different 
values for the ratio Lh/Ph- The stretching contribution to the UF state of the RNA, 
Wf, does not change when modifying the magnitude Lh/Ph, because the forces at 
which the transition occurs are quite stable under changes of Lh/Ph- However, the 
magnitude of the contribution AWj tends to notably increase as Lh/Ph becomes 
larger. 







♦ ♦♦♦♦♦♦ - 







^ 10 20 30 40 50 60 
L h /P h (contour lenght/persistence length) 

Figure 8: The different contributions to the free energy change across the transition presented as 
a function of the ratio Lh/Ph- 



6 Simulation of a pulling experiment 

As the RNA molecule unfolds-refolds in timescales much larger than the typical re- 
laxation time of handles and bead, we can consider an instantaneous relaxation for 
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these latter elements to solve the dynamical equations. This hypothesis is valid as 
long as the data is collected at frequencies smaller than the relaxational frequency of 
the bead, 10 5 Hz, that is the element with largest relaxation time (see Sec. 0J. The 
dynamics for the RNA molecule is governed by the master equation for the probability 

Pa 

d P° 1 ,L 

— = -k->p + k^pi , 
at 

^ = -k^jn + k^ Po , (24) 
at 

where fc_> and fc<_ are the unfolding and folding rates respectively. To simulate a 
pulling experiment we parallelly solve numerically the partition function of the system 
finding the mean extension and force for each element and do a numerical simulation 
of the dynamical model for the RNA 1)24(1. We implement the following algorithm: 

• We increase Xt by vAt, where v is the pulling speed, i.e the velocity at which 
the micropippete is pulled, v — Xt, and At is the iteration time, hence is 
the frequency at which data is collected 12 . 

• We compute the new (/) and (x) iteratively using the saddle point equations 
for the partition function. To these mean values we add Gaussian fluctuations 
of zero mean and variance given by (J2J). We then obtain the FEC, f(x), that 
should qualitatively reproduce the experimental one. 

• The RNA molecule is then unfolded (if it is in the folded state a = 0) or folded (if 
it is in the unfolded state a = 1) with a probability fc^(Xx)Ai and k^{X-T)At 
respectively, where At is the iteration time. These probabilities come from the 
discretization of the master equations (|24p. The unfolding and folding rates, 
k^ and fe<_, correspond to the rates for an activated process characterized by 
a barrier B(Xt) and a difference of free energy between the F and UF states 
AG(Xt) (Fig. 0(a)), 

M^t) = fc exp[-/3B(A T )] 
M*t) = fco exp[(3(-B(X T ) + AG(Xt))] , (27) 

where fco is an attempt frequency. These rates satisfy the detailed balance con- 
dition, 

J^jgj = exp[-/3AG(A T )] . (28) 

The expressions of AG(Xt) and B(Xt) are derived in appendix [B] using the 
partition function analysis. 

In the simulations presented in Secs. l6.lll6.2l we use the parameters given in Tables 
^ and [21 In Table 0] we show the values of the kinetic parameters we use, such as the 
rate of unfolding at zero force fco exp(— (3B°) and the number of opened bases n* that 
characterizes the location of the transition state (see Sec.|2J) 13 . 



12 The relation between the pulling speed v and the loading rate r (velocity at which the force increases) 
can be found using the relation between the force and displacement increments, A/ = fceff(/)AAr, as 

r = wfeeff , (25) 
where k e s is the effective stiffness of the system, computed as: 

d(f(X T )) 1 1 ! 

Ks = ^Xr- = [ k b + V x ] ' (26) 
where k x has been defined in © and kt is the stiffness of the optical trap. The F-UF transition for a small 
single domain of RNA typically occurs at forces in the range 8 — 20pN. At these forces the system verifies 
that kb is much smaller than the stiffness of the handles and the RNA molecule, kh 1 , kh 2 and k r , therefore 
we can safely assume v = r/kb- 

13 The value of n* determines the distance from the barrier to the folded conformation, xi(Xt)- With 
the assumption that n* does not depend on Xt (see footnote one can derive xi(Xt) using the WLC 
model £3 with P — Prna and L = Lrna 1 ^-, xi(Xt) = x(f) where / is the mean force acting upon the 
system when the RNA molecule is in the transition state. 
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k exp(-PB°) 


n* 


e~' M « lfr 13 
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Table 4: Parameters used to characterize the kinetics of folding-unfolding of RNA. They are 
chosen in order to reproduce the experimental kinetics results obtained with the hairpin P5ab 

In what follows we present the results of our simulations performed to analyze the 
following aspects: i) Obtaining FECs in the mixed ensemble; ii) Computation of the 
fraction of forward (reverse) trajectories that have at least one refolding. 

6.1 Force-extension curve results (FEC) 

In Figs. ITT! and II II we show the resulting FEC of our simulations for the values used 
in the experiment of Liphardt et al. shown in Tables ^ [21 an( i 01 corresponding 
to a P5ab RNA molecule and for a loading rate of r = lpN/s and of r = 50pN/s 
respectively. We do the simulation for the forward and reverse processes where Xt 
increases and decreases in time respectively. 




350 



x [nm] x [nm] 

(a) (b) 

Figure 9: Results for the FEC obtained from the simulation of a pulling experiment with r = 
lpN/s. The iteration time used in the simulation is At = 10 _2 s. In (a) we show the results of 
calculations at each time step. In (b) at we present their average over five consecutive time steps. 




100 150 200 250 
Extension (nm) 



Figure 10: Experimental FEC for p5ab obtained in JHj- The continuous line corresponds to the 
WLC curve for the handles. 
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Figure 11: Results for the FEC obtained from the simulation of a pulling experiment for a 
r = 50pN/s. The iteration time used in the simulation is At = 10 _2 s. At this pulling speed (v, 
see footnote I12JI the process is not in equilibrium and hysteresis is observed around the transition. 



As shown in Fig. at a loading rate of lpN/s different transition jumps are ob- 
served along both the forward and reverse processes, because the pulling speed (t>, see 
footnote I12J) is slow enough. In Fig. 03 (a) we represent the FEC resulting from the 
computed value of / and x at each iteration whereas in Fig. (b) we show the FEC 
obtained after averaging the results over five consecutive iterations. The amplitude of 
the fluctuations observed in Fig. 0(b) notably decreases. These values appear compat- 
ible with those found in the experimental data. Comparing these simulations results 
with the experimental FEC jlO) shown in Fig. 1101 we find a qualitative agreement, the 
shape of the curve around the transition region is qualitatively reproduced. However, 
we find some discrepancies: (i) The simulated curve is shifted in the x direction in 
comparison with the experimental one. This is because experimentally the quantity 
measured is not the absolute value of the distance x but its relative changes. Therefore 
in Fia llOl the extension represented in the a;-axis corresponds to changes in the value 
of x with respect to an initial extension of approximately lOOnm. (ii) As the force 
increases the experimental curve separates from the theoretical WLC prediction and 
therefore from the simulated results. The agreement can be improved by considering 
bigger values for the Young modulus of the handles and of the ssRNA. Furthermore, 
extending the RNA molecule model to include intermediate configurations, which de- 
pend on the number of opened bases n, one realizes that the cooperative transition 
might not be between the F (n = 0) and UF (n = N) states, but between a par- 
tially folded and a partially unfolded states. For instance, for the P5ab RNA molecule 
the cooperative folding-unfolding transition is between the state n — 3 and the state 
n = N This means that typically the first 3 base pairs open before the transition 

occurs, increasing the extension of the handles. 

Fig. II II shows the FEC corresponding to a pulling process carried out at a loading 
rate of r = 50pN/s. At this pulling speed the process is not in equilibrium and 
hysteresis effects are observed around the transition region. 

6.2 Fraction of trajectories that have at least one refolding 

We consider a system with a control parameter (generally denoted by y) that is pulled 
by changing y at certain speed v(y) — 4? . The forward (reverse) pulling process starts 
at a initial value of the control parameter yi (yj) where the RNA is in the F (UF) 
state and finishes at a final value of the control parameter yj (t/j) where the RNA is in 
the UF (F) state. We then define Np and Nr as the fractions of forward and reverse 
trajectories that have at least one refolding respectively (Fig. I12fl . 
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Figure 12: Different trajectories that have at least one refolding. The ratio between this sum 
and the total number of trajectories gives the fraction Njr (Nr) for the forward (reverse) process. 



These fractions are given by 



Np= f r< Mpvi&hja^ , ( 29 ) 

ay ay 



>i= ffMM^, (3 o) 

J Vf Jy ay dy' 

where the first integral in the right-hand side of both equations accounts for the 
probability of unfolding (folding) before a certain value of the control parameter y 
is reached and the second integral accounts for the probability of refolding once the 
RNA molecule has been unfolded (folded). The function p^ R \z, z') is the probability 
for the RNA molecule of remaining at state a until y = z' starting at y = z in the 
forward (reverse) process. The p a is solution of the master equation 

dPo {R \y,y') , , f{r), t \ , Q1 v 
j t = -k^(y)Po y '(y,y); (31) 

j t = -k^{y)Pi^'(y,y), (32) 

with initial condition p& (y,y) = 1- In appendix IC1 we prove that the fraction Np 
is equal to Nr if the perturbation protocol for the control parameter is symmetric, 
i.e. if the velocities along the forward and reverse process verify Vf(v) = —VR(y). 
In our analysis the control parameter y corresponds to the total distance Xt and 
the folding- unfolding rates are given in l|27|l . The detailed analytical expressions have 
been given (|B-5IB-6|) in the appendix[Bj However, working with these rates in order to 
do analytical computations appears quite cumbersome and it is preferable to simplify 
them. For analytical purposes we will consider effective rates where the functions B , 
AG 1 given by l|B-7(l and x\ and X2 (the distances from the F and UF states to the 
transition state along the x-axis, see Fig. |3J) are effective parameters independent of 
Xt, that we call B, AG, i\ and x%, obtaining 

fc--(/o) = fcoexp[/3(-_B + / ^i - -khx\)\ , 

M/i) - k e X p[P(-B-fix 2 + AG-^k b x 2 2 )} , (33) 

where the force f a (a = 0, 1) corresponds to the force acting upon the system at a 
given value of Xt when the RNA is in the state a 14 . In what follows we will call the 
dynamics generated by the effective rates B33JI the effective dynamics and the ones 
generated by the rates (|B-5IB-6() the non-effective dynamics. The effective rates are 
an excellent approximation to the non-effective ones in the experimental regime (see 



14 The approximation where force does not fluctuate near the transition is well justified. In fact, 

when the RNA is in a given state (folded or unfolded) the magnitude of force fluctuations is negligible (the 
r.m.s is in the range I0 -3 — fO _2 pN 2 ), so one can consider the instantaneous force equal to the mean force. 
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appendix IT))l . The relation between the forces /o and /i (for a fixed value of Xt) in 
l|3"3")l is given by 

fx = fa~ k bX m , (34) 

where x m is the distance between the F and UF states, x m = x\ + $2- Using l|34|) 
it is straightforward to see that the effective rates l(33|) satisfy the detailed balance 
condition l|28|l. We can now compute the fractions (|29I30|) as a function of the 
loading rate r. In Fig. El we show the results obtained for the fractions Np and 
Nr from the numerical computation of I|29I3U[1 using the effective dynamics (13311 . Wc 
also show the results obtained from the simulations for the fractions Np and Nr as a 
function of the loading rate r and they agree pretty well. 




r [pN/s] 



Figure 13: The fraction Np and Nr as a function of r. It is shown the results obtained from 
a 5000 realizations of the simulation of a pulling experiment and also the numerical integration 
of 129(1 or H30|) using the rates given by (|33[) characterized by the following effective parameters: 
Sln/c = 35.2fc B T, AG = 70.4fc B T, x\ = 9.75nm and x 2 = 9.35nm. 



Through the simulation we are able to compute the mean work exerted upon the 
system as a function of r: 

n 

(W(r)) = (J2fi^X T ) , (35) 
i=i 

where AXt is the increase in the total end-to-end distance in each iteration and n 
is the total number of iterations. The average is over different realizations of the 
simulation of the pulling process. The total work is the sum of the reversible work 
(i.e. the work measured in a quasi-static process for r going to zero), and the mean 
dissipated work that is also a function of r: (W(r)) — W^ v + (Wdi s (^))- 

We then consider the fraction Np for three different RNA molecules characterized 
by different parameters AG , L r , N (total number of pair bases), n* and B° lnfco and 
the results as a function of r are shown in Fig. El (a). When we plot these fractions Np 
as a function of the mean dissipated work (Wdis) exerted upon the system we see that 
the three curves corresponding to the three RNA molecules collapse to a single curve 
as it is shown in Fig. El (b) • This suggests that there is a generic dependence for the 
fraction Np as a function of (Wdis)- This dependence is not surprising as the average 
dissipated work has been already shown |18j to be a useful quantity to characterize 
the non-equilibrium regime. In particular, in the linear response regime, the average 
dissipated work depends linearly on the loading rate r, the proportionality constant 
being a function of the relaxation time of the molecule, the unfolding free energy and 
the transition force The collapse of all curves in Fig. El m a single curve is, 

however, not restricted to the linear response regime. Indeed, we have verified that in 
the regime 2ksT < (Wdis) < 5AsT, where deviations from the linear response regime 
are observable (Fig. 115(1 . there is still a good collapse in Fig. El (b) of the curves 
corresponding to the three molecules. Note that by measuring the fraction Np we 
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can obtain information about the value (Wdi s ), and knowing the total work we can 
extract the reversible work exerted upon the system. This provides an alternative way 
to derive equilibrium information from non-equilibrium experiments 15 . 
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Figure 14: (a) The fraction N F as a function of r for three different RNA molecules characterized 
by: Molecule (1) AG = 59fc B T, L r = 28.9nm, N = 24, n* = 12 and B° lnfc = 29k B T. Molecule 
(2) AG = 89k B T, L r = 40nm, N = 34, n* = 15, B°lnfc = 45fc s T. Molecule (3) AG = 39fc s T, 
L r = 16.5nm, N = 14, n* = 9 and B° lnfc = 19fc B T. (b) The fraction N F as a function of (W dis ) 
in logarithmic scale for the three RNA molecules considered in the left panel (a). 
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Figure 15: Mean dissipated work as a function of the loading rate r for Molecule 1 in (a), for 
Molecule 2 in (b) and for Molecule 3 in (c). The characteristics for the three molecules are given 
in Fig. 1141 Note that the regimes studied are far from the linear response regime as the curves 
deviate from the straight lines. The deviation from the linear response regime arises at the range 
of r where the fraction TV approaches to zero (Fig. 1141 (a)). 



7 Unfolding of domains stabilized by Mg 2+ tertiary 
contacts 

In the presence of magnesium ions {Mg 2+ ) the kinetics of the unfolding process can 
change dramatically if tertiary contacts are formed. In the experiments done in jl()j 
two different RNA molecules were studied, P5ab and P5abc, with and without Mg 2+ . 
The results obtained in ^J] show that in presence of Mg 2+ there are two different 
situations: 

• If there is no formation of tertiary contacts the folding-unfolding behavior does 
not change qualitatively. This might be consequence of the electrostatic stabiliza- 

Several methods has been proposed and tested |4UII41| 
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tion that acts homogeneously along the molecule; all base-pair hydrogen-bonds 
become more stable and the free energy landscape changes in a homogeneous 
way. This induces a slight increase of AG and B°, resulting in a value of F c 
that is a bit larger and a kinetics that is slower than in absence of Mg 2+ . Indeed, 
this is what seems to happen in the case of the P5ab RNA molecule |10| . 

• When tertiary contacts stabilized by Mg 2+ are formed the free energy landscape 
changes drastically, in particular in the vicinity of the bases that are involved 
in the formation of such tertiary contacts. Therefore the kinetics slows down 
dramatically and the unfolding-folding process changes totally, as observed with 
P5abc RNA [Kl- 
in this Section we will focus on the study of molecules that form tertiary contacts in- 
duced by Mg 2+ . Experiments on the unfolding kinetics of domains stabilized by Mg 2+ 
tertiary contacts show how intermediate states are characterized by big barriers that 
are located close to the folded state along the x-axis 16 [HAE], x\ -C x m (Fig. 03(a)). 
Consequently the height of the barrier B is quite insensitive to the force (or Xt), 
meaning that when the force exerted upon the system increases, B decreases much 
slower than the difference of free energy between both states, AG. Therefore big barri- 
ers and small values of x± imply slow unfolding processes. In complex RNA molecules 
the domains stabilized by the presence of ATg 2+ -tertiary contacts are rate-limiting for 
the unfolding of the whole molecule [331 1341 HT^l I3ri| . In these conditions, even at very 
low loading rates, the probability of refolding, once the domain is unfolded, is almost 
zero. The unfolding of RNA molecules with Mg 2+ dependent barriers at experimental 
loading rates (r « 3 — 5pN/s) becomes a 'stick-slip' process Therefore we can use 
the following transition rates 17 : 

M^t) = k e- B( - x ^/ kBT , M^t) = . (36) 

These rates have been considered by Evans and Richie in the study of bond failure 
^lEH- However, in order to get a realistic modelization using the rates l[3ti[l . the 
system requires that at the breakage force f*{Xx) (the force at which the molecule 
opens) the UF state is more stable than the F state, or AG(/*(At)) < 0. The 
breakage force changes from experiment to experiment due to the stochastic nature 
of the unfolding process. Therefore we expect that the distribution of breakage forces 
goes to zero when approaching f m , where f m verifies AG(/ m ) = 0, i.e. the value of 
the force when the RNA is in the F state at the midpoint of the transition (|14[l . i.e 
fm — /i(Ay). For such process, the kind of information that one can get from the 
analysis of the distribution of breakage forces /* is about the kinetics rather than the 
thermodynamics . 

In all the previous analysis we have considered the study of single domain RNA 
molecules. Now we want to analyze molecules that have more than one domain. 
In order to do that we extend the model developed in preceding sections to describe 
more complex RNA molecules. Here we consider how to extract kinetic information by 
analyzing the breakage force distribution for the case of a multidomain RNA molecule 
with a sequential unfolding of its domains. To this end, it is convenient to analyze 
first the case of a single RNA domain stabilized by Mg 2+ tertiary contacts. As this 
problem has been already considered by several authors we collect some of the main 
results in the appendix IT)1 

7.1 Domains with Mg 2+ - dependent barriers that unfold se- 
quentially under a loading rate 

In this section we want to investigate the applicability of the model developed for 
a single domain RNA to more complex RNA molecules such as a multidomain RNA 
molecule with a sequential unfolding of its domains. We consider a molecule composed 
by different domains under the effect of an external force, focusing on the case where 

16 Recent studies [32] show that the domains stabilized by Mg 2+ tertiary contacts are better characterized 
by kinetic models with more than one barrier. However here we will consider the simpler case of a single 
barrier per domain. 

17 Note that these rates do not verify the detailed balance condition. 
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the opening of the domains occurs in a given sequential order. There are two situations 
that favor a sequential unfolding of the domains. The first one relies on the topological 
connectivity of the molecule, that does not allow certain domains to unfold before 
others have not yet opened (Fig. El ( a ))- The second one is the blockade of the force 
induced by the most external tertiary contacts on the interior domains fFig. I16l (b)). 




(a) 

Figure 16: (a): Blockade of the force for certain domains due to the connectivity of the molecule. 
The force can not act over the domains D2 and D3 until Dl is not opened, (b): Blockade of the 
force for certain domains due to the Mg 2+ tertiary contacts. The domain D3 does not feel the 
force until the Mg 2+ tertiary contact breaks. 

For sake of clarity we will consider a sequential unfolding of a multidomain RNA 
molecule. In general the unfolding of domains is a hierarchical process not necessarily 
sequential. For instance, in left panel of Fig. El once Dl has opened, either D2 and 
D3 can be unfolded. However in our modelization we unfold sequentially the domains 
D2 and D3. The motivation to consider this simplified model is twofold. On the one 
hand, there are experimental results on the molecule L-21, a derivative of the 
Tetrahymena thermophila ribozyme, where the order of the opening of the different 
domains of the molecule studied was never observed to change. On the other hand, a 
main goal throughout this paper is to illustrate how the model for the experimental 
setup previously introduced in Sees. [21 El EJcan be generalized to include complex RNA 
molecules (and not only hairpins) rather than emphasizing details of the modeling 
of the RNA structure. With this proviso, we then model the RNA molecule as an 
unidimcnsional chain of single domains connected in series, each one represented as a 
two-states model. For a n domain system we have the F state, the UF one, and n-1 
intermediates, Ii, where i is the index of the intermediate fFig. I17(l . 
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Figure 17: Representation of the different states for a 2-domain model corresponding to a molecule 
with two domains that sequentially unfolds. The kinetics parameters of each domain are x^\ Xm 
and SW, where the super-index i = 1, 2 refers to the index of the domain. The unfolding rate and 
the breakage force for the domain i are denoted as feW and /* respectively. 

We simulate a pulling process without refolding using the effective unfolding rate 
given in (|33|l for a molecule with three domains in series. This system could represent 
the domain P4-P6 of the molecule L-21, recently investigated In these exper- 

iments, it is observed a sequential unfolding of the domains, even though there are 
different unfolding pathways because not all the intermediates are seen in each trajec- 
tory (sometimes two consecutive domains open simultaneously) . The most frequently 
observed pathway contains three transitions corresponding to the consecutive opening 
of the domains P4P6, P5 and P5abc. In Fig. ^3] we show the FEC of a 3-domain RNA 
system and in Fig. ^5] the histograms for the starting position of the rips detected. 
The results shown in Figs. El ( a ) US! (a) have been obtained by doing a numerical 
simulation of a pulling experiment using the parameters for the handles and the bead 
given in Tabled The kinetic parameters of each RNA domain are given in Fig. 1181 
In panels (b) of these figures are shown the experimental results ^Jj . 




Figure 18: Comparison of FECs between model and experiments, (a) Numerical simulations of 
the pulling process are at r = 4pN/s for three RNA domains. Simulations have been done with 
the effective model without refolding. Domains characterized by = 2.5nm flM ln(fci 1} ) = 
8.5fc B T, x^ ] = 2.5nm B^ ln(fc^ 2) ) = 8k B T, x[ 3) = 1.7nm ln(fc<, 3) ) = 8.5fc s T, where the super- 
index refers to the index of the domain. The solid lines correspond to the WLC force-extension 
curves, (b) Experimental FEC for the P4-P6 domain obtained in IJ. The solid lines correspond 
to WLC curves for the handles linked to the RNA molecule. The lower curve correspond to the 
refolding process that we do not consider here. Figure taken from jll| . 
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Figure 19: Comparison between model and experiments of the rip position distribution, (a) 
Histograms for the positions at the start of the detected rips. They correspond to the three 
transitions observed in Fig. |^| (a). The parameters used in the simulation are given in Fig. |^| 
(a), (b) Experimental histograms of rips detected in 732 unfolding curves of P4-P6 (Fig. 1181 (b)). 
Figure taken from 



For the third domain, that corresponds to the well known domain P5abc, we use 
the values of the parameters x^' and ln(fco 3 ' ) ) obtained in rjUJ. We choose the 
parameters for the other domains in order to qualitatively reproduce the experimental 
results for the unfolding trajectories in rjl] shown in the right panel of Fig. [^J The 
histograms for the positions at the start of the rips detected obtained from these 
values of parameters are different from the experimental ones (Fig. 119(1 . The main 
difference comes from the amplitude of the fluctuations of the position where each 
domain opens, that is smaller in simulations results as compared to experimental 
results 18 . Several reasons can explain this disagreement. First, there are strong 
drift effects in the machine that introduce instrumental noise. Second, no two pulled 
molecules are ever identical (disparity of the attachments, existence of more than 
one tether on the beads that can influence force measurements). A dispersion in the 
population of molecules is always another source of noise. Third, the RNA molecule 
is not just composed by a series of domains, but there are other regions (some bases) 
that do not belong to any domain. These regions can contribute differently to increase 
the length of the rips, a source of randomness for the position of the start of the rips. 
Last but not least, we cannot exclude that the kinetic model we are considering is 
too simple to explain the unfolding of these domains. It is known that complex 
RNA structures show characteristic FECs that cannot be usually interpreted in terms 
of the successive opening of native domains, because of the existence of long-lived 
intermediates including non-native helices 39 . 

The most important difference between the analysis of a single barrier (see Ap- 
pendix^ and the present study of a succession of n domains is that the force does 
not reach the domain i until the previous domain i-1 has opened. Then the domain i 
starts to be pulled only at a force larger than ff , which in our approximation is given 
by: 

ft = ft-i - hxt 1] • (37) 

where the parameters and functions with index i refer to the domain i. Let us define 
the quantity Ci as: 

a = ^§fcW(i7)- (38) 

18 Note that in experimental results the distances are given in units of nucleotides but our results are 
in units of nanometers. In the totally extended form of the polymer the conversion unit is 0.59nm per 
nucleotide. 
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The average value over different trajectories of this quantity (Ci) is a measure of 
the probability of opening the domain i just after the domain i-l has been opened. 
According to the value of (Ci) we can distinguish three different regimes: 

1. (Ci) << 1. Most of trajectories show two separated transitions (rips) for the 
opening of the domain i-l and i, because at the typical value of // there is 
a low probability of opening the domain i. It is then possible to treat the 
domain i independently of the i-l, as a single domain, using (|D-1(1 to analyze 
the distribution of breakage forces. 

2. (Ci) >> 1. The probability of opening the domain i at ff is large. Therefore 
most of the time one observes a single transition (rip) for the opening of both 
domains and the intermediate state 7j_i is hardly observable. In this case it is 
not possible to obtain information about the domain i. 

3. (Ci) « 1. This is the intermediate case between the two previous ones. We 
expect to observe trajectories with a single transition (1 rip) for the opening of 
the domains i and i-l and other ones with two separated transitions (2 rips). 
In this case, to obtain kinetic information of domain i from the analysis of the 
distribution of breakage force, /*, we need to recalculate the distribution of the 
breakage force as shown in appendix [EJ or to work with the distribution of /* 
conditioned to the fact the domain i — 1 has been opened at a force smaller than 
a given value. 

We focus on the study of the regime 3 considering two different two-domain 
molecules coupled to the system described in Sec. [21 with parameters given in Ta- 
ble 2] The system is pulled at r = 4pN/s. The first domain is the same for both 
molecules and its kinetics parameters are x^ = 2.5nm and BWh^) = 8.5k B T; 
the second domain is different for the two molecules, but both have (C2) of order of 1, 
so they are in regime 3. In order to get the kinetics parameters for the second domain, 
we will use two different techniques: 

• In appendix [El we compute the distribution of breakage forces for a domain % 
in regime 3 IJE-4J) . as a function of the kinetics parameters of domain i and 
the previous one i — 1. This technique uses the expression (|E-4j) to extract 
kinetic information for the second domain from the kinetics parameters of the 
first domain. The method consist in first building an histogram of the breakage 
forces for the second domain, using the results from all the trajectories 19 . Then 
we fit the histogram to the distribution l|E-4p 20 to get the kinetics parameters 
for the second domain. The results obtained are shown in Fig. I2"U1 

19 For the trajectories where only a single transition is observed for the opening of the first and second 
domains, it is possible to compute the breakage force for the second domain as: /| = /* — Mm ■ 
20 We truncate the series at certain k once we find convergence. 
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Figure 20: Histogram of breakage forces for the second domain of a two domain system charac- 
terized by given values of x^\ B^ ln(k^). The continuous line is the best fit to [|E-4|1 . truncating 
the series at a value k — k* where convergence is achieved. The dotted line shows the distribution 
of breakages forces for a single domain for the real values of a^ 2 ' and B^> ln(fc£ 2 ^). (a): System 



with x\ = 2.5nm and hi(k ) = 8/cbT. Series summation was truncated at k* = 2. The fit 
gives aij = 2.6±0.fnm and B^> ln(fco ) = 8.2±0.3fcsT in agreement with the correct values. The 
average value for the parameter Ci for this domain is < C2 >= 0.94. These are results obtained 
from fOOO pulls, (b): System with x\ = 1 .5nm and B' > ln(k ) = 4fc^T. Series summation was 
truncated at k* = 3. The fit gives xf> = 1.6±0.1nm and B (2 *> ln(fc^ 2) ) = 4.3±0.3fc B r in agreement 
with the correct values. The average value for the parameter Cj for this domain is < C2 >= 2.9. 
These are results obtained from f 000 pulls. 



• The second technique consist on working with the probability distribution that 
the domain i opens at a force /* conditioned to the fact that the previous domain 
opened at a force smaller than a given force < fi)- Considering 

small values of the distribution < /;) gets closer to the distribution 

of a single domain (|D- 1|> . For instance if we consider fi < f, where / is the 
minimal force at which there is no probability of unfolding the domain i at the 
given r 21 , the conditioned distribution overlaps with the distribution for a single 
domain. To compute p(f*\f*_i < fi), we do histograms of breakage forces for 
the set of trajectories that verify f*_ l < fi. Starting with a certain value of fi 
we build the histogram and do the fit to l|D-f(l to get the kinetics parameters, 
x^ and B^ \n(ko^). Then we repeat the process decreasing the value of /; until 
the parameters obtained from the fit converge to a given value; in this regime 
of values of the domain i is not influenced by the previous domain, and one 
gets the right values for the kinetics parameters. The drawback of this technique 
is that for fi too small the number of useful trajectories quickly decreases, and 
one needs many more pulls to be able to build an histogram. In the following 
figure I2T1 we show our results for p{f% \f{ < fi), for the two molecules considered 
before. 

21 The force / represents the lower limit force value below which the distribution of breakage force goes 
to zero. 
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Figure 21: Histogram of breakage forces for the second barrier of a two domain RNA molecule 
for the set of trajectories that verify /-j* < fi (the first domain opens at a force smaller than /;). 
The continuous line corresponds to the distribution of breakages forces for a single domain for the 
real values of and B^> \n(k^). (a): Same parameters as in Fig. 1201 (a,) . System characterized 
by 4 2) = 2.5nm B^Hu^) = &k B TJ t = 12pN. The fit to (fTTTfl gives if 5 = 2.5 ± O.lnm 
and \n(k^ } ) = 8.0 ± 0Ak B T. Histograms were obtained from 3000 pulls and 392 pulls verify 
/* < 12pN. (b): Same parameters as in Fig. |201 (b). System characterized by — 1.5nm, 
B& \n(kP) = 4k B T, fi = 8pN. The fit to jD^ gives xf ] = 1.5 ± O.lnm and B^ ln(fci 2) ) = 
3.7 ± 0.3fc B T. Histograms were obtained from 15000 pulls and 204 pulls verify f{ < 8pN. 



We conclude that in order to obtain kinetic information for a domain in regime 3 
both techniques are complementary The first technique has the disadvantage that it 
requires the knowledge of the kinetic parameters of the previous domain. The method 
to extract information about the second domain is to start with the analysis of the first 
domain (that is not blocked by any domain) and going forward following the sequential 
order in which domains open. On the other hand, the problem of the second technique 
is that when considering small values of the number of useful trajectories quickly 
decreases, and one needs a large number of pulls to be able to build an histogram. 
Depending on the experimental conditions one can decide which technique is the best 
to apply. 



8 Conclusions 

The recent fast development of nanotechnologies allow scientists to investigate the 
physical behavior of complex biomolecules. Of particular importance are those phys- 
ical processes in the nanoscale where the typical values of the energies involved are 
several times k B T. In such regime fluctuations and large deviations from the average 
behavior are important and deserve a careful investigation as they can contribute a lot 
to the understanding of thermal processes in small systems. RNA pulling experiments 
offer an excellent framework to address such questions as RNA molecules can be small 
enough for stochastic fluctuations be observable and measurable. 

An extremely useful technique to manipulate individual molecules are optical 
tweezers which cover a range of forces l-100pN that is relevant for many biological 
processes. A full understanding about how to extract accurate physical information 
from such experiments is therefore of great importance. The present work represents 
an attempt in that direction. At present it is not yet possible to unfold individual 
RNA molecules without attaching some polymer handles at their extremes, therefore 
all RNA pulling experiments are carried out with a system larger than the individual 
"naked" RNA molecule. This system includes the RNA molecule, the polymer handles 
and the bead in the optical trap. In order to extract accurate physical information 
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regarding the RNA molecule, a global treatment of the whole system is necessary. 

In this paper we analyzed the minimal system required to interpret the data ex- 
tracted from RNA pulling experiments. We did not include any details regarding the 
response of the machine or a realistic and accurate modelization of the structure of 
the RNA molecule. On the contrary, we have focused on those thermodynamic and 
kinetic aspects of pulling experiments by considering the transmission of the force on 
the RNA molecule induced by the the bead and handles. A key part of our treat- 
ment is a proper consideration of the ensemble that is relevant in pulling experiments. 
While the end-to-end distance (between the bead and the micropipette) and the force 
are variables that fluctuate, the total end-to end distance Xt (Fig. 01 does not. The 
thermodynamic potential in such ensemble is the key quantity that allows us to ex- 
tract accurate knowledge of the influence of these external parts (beads and handles) 
on the thermodynamic and kinetic behavior of the RNA molecule. 

In Sec. |31 we introduce the appropriate thermodynamic potential by focusing the 
analysis on single domain RNA molecules that show a highly cooperative folding- 
unfolding behavior. In Sec. [S] we analyzed the thermodynamics of the whole system 
by doing a partition function analysis that includes all parts of the setup previously 
described in Sec.0J Four are the most important results in Sec. El a) we get an explicit 
expression <|14fl for the transition force F c as well as we are able to reconstruct the 
thermodynamic force-extension curve (TFEC) from the knowledge of the parameters 
of the model, see l|20|) : b) The different contributions to the total reversible work flU} , 
coming out from the different parts of the system (bead, handles and RNA molecule) , 
have been analyzed H17I18I19|I . A comprehensive summary of these results is shown 
in Fig. 0J c) A relation between the unfolding free-energy of the molecule AGo and 
the area under the force rip W r c ip has been given in l(2*2l . d) Finally the dependences 
of the free-energy contributions to the total reversible work across the transition were 
analyzed as a function of the stiffness of the trap and the ratio between the contour 
and persistence lengths of the polymer handles (Fig. [Bland Table |2J. Taken together 
all these results establish a framework to infer thermodynamic properties of the RNA 
molecule from the experimental data. Moreover, they also allow us to understand 
under which conditions (parameters for the bead and handles) it is more reliable to 
get estimates for these properties. 

From the thermodynamics to the kinetics we verify in SecElthat the model studied 
qualitatively reproduces the results reported from experiments fFigs.l??land ll()Jl doing 
a numerical simulation of a pulling experiment. In Sec. 16. 21 we obtain some interesting 
results for other quantities that are amenable to experimental checks. In particular, 
we find a generic relation between the fraction of molecules that unfold (refold) at 
least twice during the unfolding (refolding) and the mean dissipated work. Interest- 
ingly this relation is valid beyond the linear response regime where the dissipated 
work does not increase linearly with the pulling speed. This relation could allow us to 
extract the reversible work for the unfolding process by using data extracted from non- 
equilibrium pulling experiments. This procedure is reminiscent of other techniques, 
recently applied to RNA pulling experiments |4()| . based on the Jarzynski equality or 
similar relations (for a recent review, see |41|). Moreover we have shown a symme- 
try property that relates these fractions for the forward and reverse processes. How 
general this result is in general transition state theory 021 (i- e - beyond the case of a 
cooperative two-states system) remains an interesting open question. 

In order to stress the adaptability and feasibility of our model to describe more 
complex type of molecules we consider in Sec.Hthe unfolding of a large RNA molecule 
made out of different domains that unfold sequentially. The unfolding of these do- 
mains is controlled by Mg 2+ tertiary interactions which induce large energy barriers, 
thereby a refolding event (while the molecule is pulled) is not observed at experimen- 
tal conditions. Although our study is not complete for such type of molecules (the 
assumption of a sequential unfolding may not consider other possible unfolding path- 
ways) it is instructive to see that by modifying only the model for the RNA molecule 
we are still capable of qualitatively reproducing several experimental results as shown 
in Figs. El an d El Finally we discuss possible ways to extract information about 
the kinetics of a single domain from the analysis of the breakage force distribution in 
a regime where the distribution of the breakage force for a domain depends on the 
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presence of a previous domain. 

Many aspects of RNA pulling experiments are still open, among them would be 
interesting to extend these considerations to include more complex effects induced by 
the response of the machine, test experimentally some of the results predicted in this 
work for the fraction of unfolded events and also a detailed investigation of the kinetics 
of the folding process (rather than the unfolding) in the presence of force, a process 
for which we still lack an understanding. Several of these aspects will be addressed in 
the near future. 
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A Partition function in mixed ensemble 

The partition function, Z(Xt), for the system described in Fig.Q gives the free energy 
Gx T as weu as other relevant thermodynamic properties. The state of the system is 
defined by the externally controlled variables Xt, T and P. The last two, T and P, 
are always kept at a constant value so we can ignore them throughout the paper. The 
partition function for this one-dimensional system can be written as the convolution 
of the contributions coming out from the different elements 22 : 

i°L\ pI->2 poo pL r 

Z{X T ) = C / dx hl / dx h2 / dx b / dx r Z hl (x hl )Z h2 {x h2 )Z s (x b )Z r (x r ) 
Jo Jo Jo Jo L 

xS(X T - (x hl + x h2 + x b + x r )) , (A-l) 

where Z a (x a ) is the partition function distribution of the clement a, with a = 
hi,h2,r,b. The lengths Li, L2 and L r are the contour lengths of the handles 1, 2 
and the single stranded RNA (ssRNA) respectively. The constant C is a normaliza- 
tion factor. The distribution Z a (x a ) for the element a is computed as: 

Z a (x a ) = gMe"^^ = e HJC a <»«) , (A-2) 

with j3 — j^pf- The functions E a (x a ) and G a (x a ) are the energy (or enthalpy) and 
the Gibbs free energy of the element a respectively. Both are related by G a (x a ) — 
E a {x a ) — TS a {x a ) where S a {x a ) is the entropy, S a {x a ) = ks ln(g a (x a )) and g a (x a ) 
is the density of states. We now compute the free energy G a , of each of the different 
elements at fixed value of x a : 

• Bead trapped in a potential well: As the V b (x) is the potential of mean- force for 
the bead in the trap along the reaction coordinate (see footnote [2J we can write, 

Z s (x b ) = e - pVb{Xb) . (A-3) 



Handles: We use the fact that the difference of free energy between the state 
defined with x = and the one with x — xj li is equal to the reversible work 
performed by stretching the handle from x = to x — xj li , 

GfnM = f h% dxf hi {x) = W hi (x hi ) , for i = 1,2 (A-4) 
Jo 



212, 



We restrict the configurational space to positive values of the variables x a , with a — hi,ti2,r,b. The 
reason of taking this simplification is because when considering positive values of the control parameter 
Xt the configurations with negative values of some x a have practically no weight. 
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where (x) is the thermodynamic force-extension curve (TFEC) of the handle 
i 23 . We get 

Z hi {x hi ) = e- pWh i^ . (A-5) 

• RNA: The partition function Z r can be divided in two parts, one corresponding 
to the F state (a = 0) and the other to the UF state (cr — 1). In the present 
analysis we are considering that the F state is represented by a single configu- 
ration while the UF states are represented by a continuous set of configurations 
corresponding to the difference extensions of the ssRNA (Fig. 0(b)). Therefore 
Z r is made up of two terms: a singular contribution that comes from the F 
state (cr = 0) represented by a delta function and a continuous contribution that 
comes from the UF state (cr = 1). We take the F state as the reference state with 
zero free energy. The free energy of the UF state has two terms: the free energy 
at zero force, AG , plus the corresponding loss of entropy due to the stretching: 

Z r (x r )^Z(x r ,a^O) + Z(x r ,a = l) = S(x r ) + C r e-^ AG ° +w - {x - }} , (A-6) 
where W r {x r ) is computed as in (|A-4|) 

W r {x r ) = I " dxf r {x) , (A-7) 
JO 

being f r (x) the TFEC of the ssRNA polymer. The probability P(a) for the RNA 
molecule to be in the state a is given by -P(cr) oc f^ r dx r Z{x ri cr). To compute 
C r we use that the RNA molecule at zero force satisfies 



^ = -^M^) (A- 



and substituting (IA-6|) we obtain, 



C r = —j . (A-9) 

J Q Lr dxe-WrW 



Adding the different contributions we get: 

pLi poo pL r 

Z[Xt) = CI dxhx I dxh 2 / dxb / dx r 
Jo Jo Jo Jo 



-P(W hl +W h2 +V b ) 



x [S(x r ) + C r e-^ AG " +W ^]6(X T - (x hl + x h2 + x b + x r ))] . (A-10) 

We now separate (| A- 10|) in two contributions coming from the F and the UF states. 
By using the integral representation of the delta function, 

1 [°° 

8{x) = — / exp(iAx)dA (A-ll) 



2tt _ 
we get 

Z{X T ) = Z (X T ) + Z X {X T ) with (A-12) 

Z Q (X T ) = ^ rfA e (lAXr+9o(A)) and ZAX T ) = / dAe (lAXr+9l(A)) ,(A-13) 

27T /_„ Z7T /_„ 
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Which variables are controlled is a relevant choice in single molecule pulling experiments. In contrast 
with macroscopic systems, where all the ensembles are equivalent |15|. FECs depend on the particular 
ensemble considered. Eq. has been defined for the isometric ensemble. The isometric TFEC is 

the thermodynamic curve corresponding to a system in the ensemble where the end-to-end distance x is 
held fixed, and is given by the mean force as a function of x, (f)(x). While the isotensial TFEC is the 
TFEC resulting of working in the force ensemble, (x)(f). In general both TFEC differ |15) . but in this 
analysis we consider that the handles and the RNA molecule are long and flexible enough to have an 
identical isometric and isotensional TFEC that we call f a (x a ) with a — h\,h2,r. This allow us to use the 
extrapolation expression JJJ (or the one given in for the function f a (x a ) when using the WLC model 
to describe the polymer behavior. 
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where the functions go and g\ are given by 



.9o = log 



dx 



In 



dxh 



dxb\e 



-P(W hl +W h2 +V b ) -i\(x hl +x h . 2 +x b )l 



(A-14) 



9i = log 



L\ poo pL 1 

dx hl / dx h2 / dx b 
o Jo Jo Jo 

e -i\(x hl +x h , 2 +x b +x r )l 



dx r [C r e~^ w ^ +Wh - +H+AG°+W,) 

(A-15) 



Eqs. \ A- 1.3(1 for Zq and Zi are integrals respect to A of an exponential with an argument 
that is extensive with the size of the system 24 . Therefore if the system is big enough, 
the saddle point approximation is valid and becomes exact in the thermodynamic limit. 
As a check we have verified that the results from the saddle point approximation and 
the exact numerical integration of the partition function are in pretty good agreement 
for the system with parameters given in Tables ^ an d [3 Applying the saddle point 
technique, one is led to extremize the arguments of the exponentials with respect to 
all the variables of integration. In this way we obtain: 



A CT with cr = 0,l and a = hi, hi, r, b , 



(A-16) 



where x^ corresponds to the value of the variable x a when the RNA molecule is in 
the state a that extremizes the argument of the exponential. We have two branches 
corresponding to the situations where the RNA is folded (a = 0) and where the RNA 
is unfolded (a = 1). We use the super-index a to denote each branch. Eq. I|A-16|) 
tells that the integration variable A plays the role of the thermodynamic force, so the 
Act corresponds to the mean force acting upon the system for the branch a and for a 
fixed value of Xt called (f) a . Eq. jA~16)l can be written as 



/?(*2) 



/ft 1 (4 1 ) = 4(4j = /r(^) = (/>l 



(A-17) 



where the force f° 



, dW°(x) y 

dx 



is the mean force acting upon the element a at fixed 



x a = x a a for the branch a. In Fig. 1221 (a) we show the two branches (f) a as a function 
of Xt for a system with parameters given in Tables Q and HI The transition from the 
F-UF states is the jump from one branch to the other. 



24 By size we mean the length of the handles as well as the length or molecular weight of the RNA 
molecule. In general to apply the saddle point approximation we require that the energies of the different 
elements of the system (bead, handles and molecule) are several times ksT. 
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(a) 



300 400 
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Figure 22: We consider a system with the parameters given in Tables^and El From the partition 
function analysis we compute: (a) The two branches (/) CT , corresponding to the thermodynamic 
forces acting upon the system for a given a RNA state as a function of Xt- (b) The free energy 
Gx T an d the free energy of each branch c, G a , as a function of Xt- 



Hence we have obtained that the values of the arguments for which the contribution 
to the partition function is maximum corresponds to the equilibrium values at a given 
Xt: 

Z(X t ) = Z (X T ) + Z X {X T ) , (A-18) 



Z (Xt) «exp -l3(W hl ({x hl ) ) + W h2 ((x h2 ) ) + V b ({x b ) )) 



(A-19) 



Zx(X T ) « exp[-/3(W hl ((a: hl )i) + W^xQx) + V b ((x b )i) + AG + W r {{x r )x)) . 

(A-20) 

where we have neglect the subdominant contributions. The {x a ) a correspond to the 
mean value of x a for the branch a and for a fixed value of Xt- In Fig. [221(b) we show 
the results for the free energy of the system with parameters given in Tables ^and [21 
as a function of Xt, 

G Xt = -k B Thx{Z{X T )) , (A-21) 



and also the free energies of the system for each branch a, 

G a = -k B T\n(Z a {X T )) • 



(A-22) 



The free energy of the system Gx T changes from one branch to the other at X T , when 
both states are equal probable, i.e Go = G\. 



B Computation of the folding and unfolding rates 
in the mixed ensemble 

We model the kinetics of the folding-unfolding of RNA as a Kramers activated process 
characterized by the following transitions rates: 

M^t) = fc exp[-/3B(A T )] 
fc_(A T ) = fc exp[(3(-B(X T ) + AG(X T ))} , (B-l) 

where ko is an attempt frequency that depends on the shape of the free energy land- 
scape, on the molecular damping and on the natural frequency of the hydrogen bond 
oscillations ^21- The functions AG(Xt) and B(Xt) are the difference of free energy 
between the F and UF states and the height of the kinetic barrier located between 
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them (Fig. 13(a)) 25 . Using the results obtained from the partition function analysis 
we can write AG(Xt) as: 

AG(Xt) = -k B T\n Zl( -* T \ = AG" + W r {(x r ) x ) - (f) x m + \k b x 2 m + AW h ,(B-2) 

where we used (|9I10|) for the expressions of Zq and Z\ and we used the parameter x m 
defined as the distance between the two states, x m — (x)i — {x)q. The functions W r 
and AW h are given by JHJ and 

The height of the barrier is given by the difference of free energy between the F state 
an the transition state that we will denote as a — t (averages taken when the molecule 
is in its transition state will be denoted by (...)t). The transition state is located at 
the point where the free energy landscape of the system depicted in Fig.Q]is maximum 
(Fig. El(a)), and we define it as the RNA state where the first n* bases are opened and 
the latter N — n* are closed, ./V being the total number of bases that form the RNA 
molecule. Therefore the function B(Xt) is computed as the free-energy difference 
between the folded state F and the transition state, which are separated by a distance 
x\ = (x) t — (x)q. This gives 

B(X T ) = B° + W r ((x r ) t ) - (f) x 1 + \/2k b x\ + AWl . (B-3) 

The function W r is given by and AW^ is the change in free energy of the handles 
when the RNA molecule jumps from the F state to the transition state computed as: 

AW t h = W hl ((x hl ) t ) + W h2 ((x h2 ) t ) -W hl ((x hl ) Q ) -W h2 ((x h2 ) ) . (B-4) 

Then the rates fc_» and fc<_ associated to the activated process can be written as: 

M*t) - fc exp[/3(-S 1 + (f) oXl - l/2k b x\)] , (B-5) 

M^t) = fco exp[/3(-B 1 + AG 1 - {f) ± x 2 - l/2k b x 2 2 )] , (B-6) 
with B 1 = B° + W r ((x r ) t ) + AWl , AG 1 = AG + W r {{x r )i) + AW h , (B-7) 

where we used (|B-1|) . (|B-2|) and (|B-3|I . The expression for the rates <|B-5IB-6t are 

equivalent to the ones obtained by Bell jHJ but in the mixed ensemble. Note that the 
two rates fe_ k (Xx), fc<_(Xr) satisfy the detailed balance condition l)28[l. 

C Demonstration of the equivalence between the A^^ 
and Nb 



Taking the expressions for the fractions Np and Nr given by I|29l3t)|l and integrating 
the left integral we get: 

N F = 1- pi( yi ,y f ) + df>F °^ y) p[(y,y f )dy , 



Vi 



dy 



Vl dpf(yf, y) 
dy 



N R = 1- p?(y f , + / ^^/"' pZiy, yi)dy , (C-1) 
Jvt 



where y denotes a generic control parameter. Then using the equation for the evo- 
lutions of the probabilities p a given by (|31I32|I and for a symmetric perturbation 



protocol, v F {y) = -gf = -v R (y) - w 



'Hi 



, we obtain the following relation: 

R 



it (//'-//> = -xp - /" ^yW i =expf- [ V K ^'}yP d y "]=p?(y,y>)(C-2) 



v' v F (y") 1 I J y v R (y") 



25 Note that the physical meaning of AG(Xt) is completely different from AGx T (see 115> ). The latter 
corresponds to the free energy difference of the global system between two different values of Xt- 
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We consider the expression for Nr given by IjC-lJl and we integrate by parts, 



N R = i- Pi(yf,y t ) - Poiyj, yi) + pi(yf, yi) 



y.f 



dy 



P?(y f ,y)dy. (C-3) 



Using the relation between the probabilities p a for the forward and reverse process 
(jC-211 . we obtain 



N R = l-pg{vi,yf) + 



vs dp%{y t ,y) 



pf(y,yf)dy = N f 



(C-4) 



D Single domain RNA as a stick-slip process 

To address the case of a multidomain molecule it is useful to focus first on the simpler 
case of a single domain molecule. We consider an unfolding process without refolding 
(fc^_ = 0) characterized by an effective unfolding force-dependent rate k_(f) J23J. 
The distribution of breakage forces is given by ^2 



-[fc^(/*)-k^(0)] 



(D-l) 



Note that in order to the no refolding condition to be realistic there must be a limit 
force f m , below which the distribution of breakage forces goes to zero. This lower limit 
f m arises because we are considering that there is a vanishing probability of jumping 
if the UF state is not thermodynamically stable, AG(f m ) — 0. From this distribution 
one can compute the mean value and the variance of the breakage force 



(/*) = / dfP(f)f 



Xi 



-e a Ei{-a)} 



(D-2) 



where a = ^|^/c^(0) and the function Ei is the special elliptic function. By doing 
an expansion in the parameter a, that is much smaller than one (otherwise there is a 
finite probability of refolding) , we obtain: 



</*> 



k B T\ 



In 



^k B Tk e-^ B -^ 2k ^ 
where 7 is the Euler's constant, and 

,k B T 



0(a) 



4 



(/"> - (tr 



(^)V/6] 



0(a) 



(D-3) 



(D-4) 



Therefore by studying either the distribution of /* at fixed r or the mean value or 
the variance of such distribution as a function of r one can obtain information about 
the kinetic parameters doing a fit to i|D-l|) . (|D-3|) or i|D-4|) respectively 26 . 

In Fig. 123 we plot (/*) and a"j, as a function of r obtained by pulling the system 
described in Sec. with parameters given by Table ^ The kinetic parameters that 
characterize the RNA molecule that we consider here are given in Table El 



k exp(-pB°) 


n* 


e~ y w lO" 4 


2 



Table 5: Parameters that characterize the unfolding kinetics of the RNA hairpin. 



We perform two kinds of simulation both using the condition of no refolding l|36|l , 
but with the dynamics generated by different unfolding rates, the non-effective rates 
l|B-5(l and the effective rates (|33|l respectively: 

26 The system under consideration verifies that the transition occurs close to the situation where the kb 
is much smaller than the other stiffness values, ktn, kh 2 an d k r . Beyond this regime, one should take into 
account the variability of r with the force. And to get a more accurate result one should fit the data to 
the distribution of the breakage forces instead of the mean or variance of the breakage force as discussed 
in [22]. 
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• Non-effective rates: We consider the explicit dependence on Xt of the barrier 
B(Xt) that governs the unfolding kinetics, see l|B-3j) . 

• Effective rates: We use the unfolding effective rate ll-iiill in order to generate the 
unfolding dynamics where we neglect the dependence on Xt of x% and B 1 as 
given by 1B-7J1 . For kinetics processes with barriers quite insensitive to the force 
(or Xt) this seems to be a reasonable approximation. The effective model is also 
the one we use to do the analytical computations. 

The comparison between both simulations allows us to see how big are the differ- 
ence between both models, and how far the analytical results are from the non-effective 
model. In Fig. |2S1 (b) we show a"j, for both kinds of simulations for a broad range of 



values of 



The non-effective simulation gives fluctuations aj. 



that decrease when 



increases, instead of being constant as the effective model predicts. This effect comes 
from the dependence of x\ on / (or Xt)- On the other hand we can see that when 
r approaches zero, fluctuations disappear, because the domain is always opened at 
zero force 27 . However, we should note that in the latter regime (r going to zero) the 
no-refolding approximation becomes invalid, because the distribution P(f*) do not 
vanish for /* < f m . Nevertheless we see that for the range of interest r « 1 — 50pN/s 
both simulations agree pretty well either for the ex 2 , as for the (/*). We can conclude 
that the effective dynamics reproduces well the non-effective one. 
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Figure 23: Comparison between the effective and non-effective rates. Results for 3000 pulling 
trajectories for a domain characterized by the kinetics parameters given in Table |SJ For the 
effective kinetics parameters we use x\ = 2.5nm and £?ln(fc ) == 8.5/esT. In (a) it is shown the 
mean breakage force as a function of the pulling rate. The straight line is the best fit to a function 
y = ahx(x) + b, obtaining x\ = 2.47 ± 0.03nm, Z?ln(fc ) = 8.3 ± 0.2fcsT. In (b) we represent the 
variance in breakage force as a function of the pulling rate. The straight line is the best fit to a 
constant y — C for the data with r > lpN/s, obtaining x\ = 2.45 ± 0.06nm. 



Fitting the data obtained from the simulation for (f *(r)) to <|D-3|) . and for er 2 , to 
l|D-4(l . we get accurate results for both parameters, x\ and B\n(k ). 



E Computation of the distribution of probability of 
breakage force in regime 3 



According to (|D-lJl the distribution of breakage forces for the ith domain conditioned 
to a given value of ff (with // = f*_ x - k b xt' 1} ), p{f*\lt), is: 



27 In the simulation there is no restriction for the breakage force, hence /* can be smaller than f m . When 
r goes to zero the breakage force too, because the molecule always opens if we wait long enough (and does 
not close anymore as = 0). 
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P(/*I//) = OhU* - f!) , (E-l) 

where the parameters and functions with index i refer to the domain i. The Oh is the 
Heaviside function, Oh{x) = 1 only if x > otherwise Oh(x) — 0. Assuming it -1 ' 
as a constant parameter we derive the breakage force distribution P(f*) averaging 
(|E- 1 f) over the distribution of the breakage forces of the previous domain i-1, P{f*_i)- 
To get P(/*_ 1 ) one has to average over the distribution of the breakage forces of the 
domain i-1, and so on. This leads to the following recurrence formula: 

/CO 
p(mfi) n <Vkp(fi\fkWiP(fi) > (e-2) 

-°° k=2 

where C is a normalization factor. We will consider the case where the distribution of 
breakage force for the domain i-1 is not modified by the previous one, either because it 
is the first domain, or because the typical value of f*_ x is higher than all the previous 
ones with k < i — 1. Then the distribution of f*_ 1 is computed as in the case of a 
single barrier HD-ll) and (|E-2I) reduces to, 

/>oo 

P(f*) = C / P {mf!)P{fU)dfU . (E-3) 



The integral in i|E-3() is not analytically solvable. Then we expand in power series 
the exponential exp — [k^ (//)] in (|E-1I) and substituting in (|E-3|) we obtain: 



P{ft) - E f ] , (E-4) 



where 



k B T 



A\ = \ QP(/ti)C-i 

— OO 



[(feT/raf" 1 ')^ 1 ^)]^"^^ 



x 7 (^i W M ( ^ 1) + 1, — Tr7T^(/; + • (E-5) 

The j(x, y) is the incomplete gamma function. For the regime 3 the series can be 
truncated, because the moments of Ci, M, are not big and the series fastly converges. 
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